Dynamical vortex phases in a Bose-Einstein condensate driven by a rotating optical 
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We present simulation results of the vortex dynamics in a trapped Bose-Einstein condensate in 
the presence of a rotating optical lattice. Changing the potential amplitude and the relative rotation 
frequency between the condensate and the optical lattice, we find a rich variety of dynamical phases 
of vortices. The onset of these different phases is described by the force balance of a driving force, 
a pinning force and vortex-vortex interactions. In particular, when the optical lattice rotates faster 
than the condensate, an incommensurate effect leads to a vortex liquid phase supported by the 
competition between the driving force and the dissipation. 
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Driven vortex lattices interacting with periodic pin- 
ning sites have attracted considerable interests, because 
they exhibit novel phase transitions between dynami- 
cal vortex states as a function of driving force. This 
problem has been studied in a wide variety of con- 
densed matter systems including type-II superconductors 
and Josephson-j unction arrays, described by the Frenkel- 
Kontrova-type models of friction [1]. Recently, alkali- 
atomic Bose-Einstein condensates (BECs) in periodic po- 
tentials have been proposed to be another possible, su- 
perior system to study the vortex state caused by the 
interplay between pinning effects and vortex- vortex in- 
teractions 0, 01 51, IE @] . In the condensed matter systems 
such as superconductors, the pinning effects incidental to 
unavoidable impurities play an essential role in the vortex 
dynamics. Since the BEC system is free from impurities, 
the intrinsic phenomena associated with quantized vor- 
tices can be studied through direct visualization [7|]. In 
addition, a laser-induced optical lattice (OL) constructs 
an ideal periodic potential, in which the microscopic pin- 
ning parameters, such as depth, periodicity, and symme- 
try of the lattice can be carefully controlled 0. 

In this Letter, we study the dynamic properties of vor- 
tices in a trapped BEC driven by a rotating OL. A ro- 
tating OL has already been realized experimentally by 
using a laser beam passing through a rotating mask [6[. 
In the previous studies [J 0, 0, E[ , the authors consid- 
ered the equilibrium properties of a rotating BEC in the 
presence of a corotating OL with the same rotation fre- 
quency as the BEC. They obtained the phase diagram 
of the vortex lattice structure as a function of the pin- 
ning parameters through variational Q or numerical 0] 
analysis. For a strongly interacting limit, this system 
can be mapped into the Bose-Hubbard Hamiltonian in 
the rotating frame 0, [H[ , exhibiting a series of quantum 
phase transition Q . Tung et al. reported the experimen- 
tal demonstration of the vortex pinning and the evidence 
of the structural phase transition [6[ . Beside these inter- 
esting studies, we can also consider a problem in which 



the rotation frequency of the OL is different from that of 
the condensate. Then, the moving optical lattice yields 
a rich variety of nonequilibrium vortex phases that are 
characterized by distinct vortex dynamics. The onset 
of these phases can be understood by the simple argu- 
ment of a balance between the driving force, the pin- 
ning force, and the force coming from vortex- vortex in- 
teractions. In particular, when the optical lattice rotates 
faster than the condensate, the competition between the 
driving force and the dissipation causes a vortex liquid- 
like phase. This is an example of dissipative structures 
sustained by energy input and output [9[ , which offers an 
avenue to study them in atomic-gas BECs. 

This system allows us to run ab-initio simulations 
based on the Gross-Pitaevskii (GP) equation fioj ] . 
whereas the other system has relied on phenomenolog- 
ical theory [1]. Suppose that the condensate is created 
in a rotating harmonic trap with an angular velocity Qz. 
Here, we confine ourselves to the pancake geometry by 
assuming a tight axial confinement of a harmonic trap 
and high rotation rates which reduce the radial trap fre- 
quency co±. We then give the Gaussian for the frozen 
axial component of the condensate wave function, obtain- 
ing a quasi-two-dimensional GP equation with a renor- 
malized coupling constant. In a frame rotating with Q, 
the transverse wave function ip(x,y,t) (normalized to be 
unity) obeys 
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we choose our units for time, length, and energy as 
1/uj±, a± = ^jTi/2muj^, and hu±, respectively. Here, 
r 2 = x 2 + y 2 , m the atomic mass, \i the chemical po- 
tential, and L z the z component of the angular mo- 
mentum operator. The coupling constant is given by 
U = 4VirXNa/a± with the particle number iV, the s- 
wave scattering length a, and the aspect ratio A = uo z /uo±_ 
of the trapping potential. The pinning effect is caused 
by the OL potential Vol = ^o[sin 2 (/cX) + sin 2 (kY)} with 



2 



(a) 



(b) 8co =-0.1 



Sco=0 



8co=0.1 



* * $ 



(c) 



o 4 







1 


i i 




£7=1000 
^2 = 0.85" 




fully pinned 
vortex lattice 




7= 0.01 


inner-pinned \ 
outer-depinned\ 
Vortex lattice * 










\ pre -melting 
vortex lattice. 




vortex 
liquid 


vortex lattice 
sliding , 


i - - i - - ' 


v.ortex lattipe sliding 



-0.3 -0.2 -0.1 







0.1 0.2 0.3 



FIG. 1: (a) The profile of the OL for k = tt/3 at t = 0. Bright 
(Dark) regions represent the potential maxima (minima), (b) 
Condensate density for each Suj at a certain time in the quasi- 
stationary stage. The parameter values are U = 1000 and 
Q = 0.85. The middle panel (Suj — 0) is the initial state of 
the simulations. The presented box dimensions along x and y 
are —11.8 to +11.8 (The simulations were done in the region 
— 14.0 to +14.0). (c) Suj (relative rotation frequency) - Vb 
(potential amplitude) phase diagram of the dynamical vortex 
states. The phase boundary is determined by calculating the 
structure factor and the angular momentum per atom, as de- 
scribed in the text. The dashed boundary for the depinning 
of all vortices is obtained by solving Eq. J4}. 



square symmetry. We consider the situation in which the 
OL rotates with a rotation frequency uj different from Q 
of the condensate. Then, Vol has an explicit time de- 
pendence through 



cos Sujt sin Sujt 
— sin Sujt cos Sujt 



(2) 



Here, Suj = uj — Q is the relative rotation frequency be- 
tween the condensate and the OL; the OL rotates faster 
(slower) than the condensate for Suj > (Suj < 0). Fi- 
nally, we introduce the phenomenological parameter 7 in 
the left side of Eq. (pQ), which represents the damping 
of the condensate motion caused by the noncondensed 
atoms rotating with ft [11]. The damping is necessary to 
obtain the phase diagram of the quasi-stationary states. 
We use the small value 7 = 0.01 which is small enough 
to ensure the underdamped motion of the vortices; the 
overdamped motion is unlikely to occur in an ultracold 
dilute atomic BEC. 

Our simulations start from the initial state in which the 
BEC is equilibrated under rotation Q and a corotating 
OL with the same rotation frequency uj = Q. Although a 
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FIG. 2: (Color online) Time development of the structure 
factor 5(ksQ, £) for Vb=2.0 and for different values of Sou. The 
average of the obtained data is given by the solid curve. The 
inset show the corresponding vortex trajectories within the 
Thomas-Fermi radius of the initial state in a frame rotating 
with the OL. 



rich variety of equilibrium vortex structures emerge de- 
pending on the system parameters [3[, we focus on the 
initial state in which vortices are commensurate with the 
pinning sites. For example, given Q = 0.85, the wave 
number k = tt/3 for the OL [see Fig. QJa)] and the cou- 
pling constant U = 1000, we can obtain such an initial 
state through the imaginary time propagation of Eq. (pQ), 
as shown in the middle of Fig. QJb). The nucleated vor- 
tices are pinned at the potential maxima of the OL and 
form a square lattice. We do not consider the very small 
values of Vq either, because the competition between tri- 
angular vortex lattice and the square OL potential occurs 
0. Next, we tune the rotation frequency of the OL by 
Suj adiabatically so as not to excite the collective motion 
[I2I ]. We then calculate the time evolution from this ini- 
tial state for a given set of parameters (5w, Vb), using the 
Crank-Nicolson implicit scheme. 

After a sufficiently long time, the condensate motion 
becomes quasi-stationary because of the damping pa- 
rameter 7, with the remaining durable vortex dynamics 
driven by the OL. These dynamics exhibit nontrivial be- 
havior depending on the values of (5w, Vb), summarized 
in the Suj - Vo phase diagram in Fig. DJc). To characterize 
clearly the dynamics and the underlying structures, we 
calculate the vortex trajectory in the frame rotating with 
the OL and the corresponding structure factor defined as 
5(k, t) = jjr- ^\ exp[zk • Tj(t)] [3|, where j labels individ- 
ual vortices, Yj(t) is the position of the jth vortex and 
N c is the total number of vortices. If vortices are pinned 
by the OL and form a perfect square lattice, 5(k, t) be- 
comes unity if we choose k = k$Q = 2fcy, which is one 
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FIG. 3: Time- averaged values of the structure factor (S'(ksQ)) 
(a) and the angular momentum per atom (l z ) (b) as a function 
of 5uo for several values of Vb. In (a), the result for 7 = 0.05 
and Vb = 1.0 is also shown. In (b), the plot of the rigid-body 
estimation [13] l z — C(Q + Sou) /y/l — (fi + Sou) is shown by 
the dotted curve; C — 6.7 is the fitting parameter. 

of the two fundamental reciprocal lattices for the square 
lattice. Figure [2] shows the time evolution of 5(ksQ, t) as 
well as the vortex trajectories. Here, we pick out the vor- 
tices inside the Thomas-Fermi radius of the initial state 
to obtain Fig. O The time-averaged value of 5(ksQ,t) 
and the angular momentum per atom l z are the useful 
quantities to characterize each distinct phase, which are 
shown in Fig. [3j 

There is a remarkable difference in the dynamics be- 
tween the positive side of Sou and the negative one. For 
Sou < 0, the rotation of the OL is slower than Q. If 
\6u\ is small, the initially pinned vortices stay in their 
original pinning sites, being dragged by the rotating OL. 
This behavior is shown in the inset of Fig. E^a); 5(ksQ) 
stays about 0.8 in the stationary stage. Then, the con- 
densate radius and the angular momentum of the con- 
densate decreases [see Fig. D^b) and Fig. [3fb)] because 
its rotation must be reduced by the dragged vortices. An 
increase in \5w\ causes depinning of vortices from the con- 
densate periphery as shown in Fig. Efb), which leads to 
the decrease of the time- averaged value (5(ksQ)). With 
further increase in \5u\, the vortices cannot catch up to 
the motion of the OL and start to make a sliding mo- 
tion. The onset of the sliding is seen from the devia- 
tion of (l z ) from the rigid-body estimation. When the 
sliding motion occurs, all vortices rest in the rotating 
frame of Q. We identify the vortex phase as follows: 
(i) "fully pinned vortex lattice" which is determined by 
(S'(ksQ)) > 0.7, (ii) "inner-pinned outer-depinned vortex 
lattice" by 0.25 < (#(k SQ )) < 0.7 and by (l z ) following 
the rigid-body estimation, and (iii) "vortex lattice slid- 
ing" by (5(ksQ)) < 0.25 and by (l z ) deviating from the 
rigid-body estimation. 

For 5uj > 0, if vortices are pinned, the angular mo- 
mentum of the condensate must increase because the OL 
rotates faster than Q. This demands additional vortex 
nucleation to ensure the fast rotation. The nucleated 
vortices give the incommensurate effect for the vortex 
lattice. These vortices wander interstitially, exchanging 



their locations with the pinned vortices, where the vor- 
tex trajectories show that the vortices move to nearest- 
neighbor sites. For small 5(j, the vortex configuration 
keeps the original square symmetry, (5(ksQ)) staying 
about 0.85. However, a slight increase in 5uj rapidly de- 
creases (5(ksQ)), which indicates depinning of some or 
all commensurate vortices. We find that the interstitial 
vortices induce melting of the vortex lattice from the in- 
side as shown in Fig. [2fc). If Vb is relatively small, (l z ) 
in the quasi-stationary stage is near the initial value and 
there are few interstitial vortices. The vortices then make 
a sliding motion. For large Vb on the other hand, the sta- 
tionary value of (l z ) is larger than the initial value and 
the incommensurate effect becomes important. Then, 
as shown in Fig. [2jd), the vortex motion is fully disor- 
dered by the nucleated additional vortices which prevent 
the other vortices from pinning, exhibiting a liquid-like 
behavior. We characterize the phase (i) "fully pinned 
vortex lattice" by (S'(ksQ)) > 0.85, (iv) "pre- melting 
vortex lattice" by 0.25 < (S(k SQ )) < 0.85, (v) "vor- 
tex liquid" by (S(k S Q)> < 0.25 and (l z ) - l z (t = 0) > 2, 
and (vi) "vortex lattice sliding" by (S'(ksQ)) < 0.25 and 
(l z )-lz(t = 0)<2. 

The onset of these different phases can be described 
using force balance arguments. The force acting on a 
vortex is given by F = + F p + F vv , where F^ is the 
driving force, F p the pinning force, and F vv the force 
from the other vortices. The commensurate vortex will 
remain pinned as long as the inequality \F p \ > | F^ | + 1 F vv | 
is satisfied. The pinning potential /7 pin for the vortex 
displacement r = (x,y) was calculated by Reijnders and 
Duine [2] for the case without a harmonic trap. Using 
their results, we obtain the magnitude of the pinning 
force F p = — V/7 p i n as 

|F P | = ^kV Q(K)^m 2 (2kx) +sin 2 (2%), (3) 

where a z is the size along the z-axis, £ = (87rap) _1//2 
the healing length with the condensate density p, and 
Q( k = JT drJ (2kr)/r + Ji (2/c£)/2/c£ with the Ith- 
order Bessel function J\ of the first kind. This analytic 
form is valid when the vortex core is much smaller than 
an optical lattice period and the strength of the potential 
is sufficiently weak. 

For Su < 0, because the incommensurate effect is 
not important, the square symmetry of the pinning site 
makes the F vv irrelevant for the vortex dynamics. The 
driving force is given by the sum of the Magnus force 
Fm and the mutual friction force F^; the mutual fric- 
tion originates from the interaction between the vortex 
core and normal-fluid flow. If a vortex is pinned by the 
rotating OL, the vortex experiences the background flow 
with — 5uo0 in the rotating frame of w. Then, the Mag- 
nus force is written as Fm = a z mpKSujrz x 6 with the 
quantum circulation k = h/m. Since the Magnus force 
becomes larger as r increases, vortices begin to depin 
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from the condensate periphery. At near zero tempera- 
ture, the thermal atoms and the resulting mutual fric- 
tion force can be negligible. Hence, the sliding motion 
can occur when |Fp| for vortices near the trap center 
is less than |F^| on the same position. We focus on a 
vortex located at the pinning site (x,y) = (7r/2fc, ir/2k). 
By comparing the Fm at r = n/V2k and the maximum 
of |F P | at (x,y) = (7r/4fc, ir/Ak), we can obtain the criti- 
cal 5(j for the depinning of all vortices. To estimate the 
chemical potential /i and the healing length £, we use the 
centrifugal-force modified Thomas-Fermi approximation 
that neglect the density variation caused by the vortex 
cores. Using the central density p = m/i/Airf^a, the crit- 
ical 8u) can be obtained by solving the equation 

h"VoQ(k^)) = J W-^-^ Su,. (4) 

7T V Z7T 

The result is shown by the dashed line in Fig. [He), which 
is in reasonable agreement with the numerically obtained 
boundary for the sliding motion. 

For Sou > 0, the argument becomes more complex be- 
cause the interactions caused by the interstitial vortices 
cannot be negligible; the solution of Eq. (|4|) overestimate 
the boundary of the vortex-liquid phase because of the 
neglect of F vv . Also, the Thomas- Fermi approximation 
breaks down for Q + Sou = uo c± 1 (5u ~ 0.15). We need 
further detailed consideration to obtain the onset of the 
vortex lattice melting caused by the incommensurate ef- 
fect. 

The melting of vortex lattices is one of the most im- 
portant topics in atomic-gas BECs. The melting was pre- 
dicted to occur in the limit of fast rotation and low con- 
densate particle number through thermal [14] or quan- 
tum fluctuations [15]. The melting mechanism shown 
in this paper is different from those studies; the vortex 
liquid phase is supported by the driving force and the 
dissipation. This is one of the examples of dissipative 
structures sustained by the balance between energy in- 
put and output [9[. 

Using different values of the dissipation parameter 7, 
we find no qualitative change in the results if the value 
of 7 is sufficiently small. However, with increasing 7, we 
find that the region of the fully-pinned lattice is reduced 
from the 7 = 0.01 results; an example of (S'(ksQ)) with 
7 = 0.05 and Vo = 1.0 is shown in Fig. 02a). This is due 
to the increase of the driving force associated with the 
mutual friction force F^; the damping parameter 7 is 
related to the interaction between the condensate atoms 
and the thermal atoms [Tlj . 

In conclusion, we study the dynamical vortex phases 
of a trapped BEC in the presence of a rotating OL with 
a different rotation frequency from the condensate. We 
find several nonequilibrium dissipative structures, where 
the vortices exhibit complex dynamics due to the com- 
petition of the driving force, the pinning force and the 



vortex- vortex interaction. When the rotation rate of the 
OL becomes higher than that of the condensate, the addi- 
tional nucleated vortices destroy the commensurate vor- 
tex lattice, leading to the vortex liquid-like phase. 
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